import numpy as np
import matplotlib.pyplot as plt
import os
import datetime
from refract import *

# Separate into months.
months = [ 'JAN', 'FEB', 'MAR', 'APR', 'MAY', 'JUN',
           'JUL', 'AUG', 'SEP', 'OCT', 'NOV', 'DEC' ]
ndays = [ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 ]
frequencies = np.array([ 2100.0, 5500.0, 9000.0, 17000.0, 22000.0, 33000.0, 43000.0, 83000.0, 100000.0 ])
hfrequencies = frequencies * 1000000.0
linecols = [ 'black', 'green', 'red', 'blue', 'cyan',
             'magenta', 'yellow', 'purple', 'orange' ]

startYear = 2010
endYear = 2014

loadNeeded = False
for m in xrange(0, len(months)):
    fname = "mediantemperature_%s.txt" % months[m]
    if (not os.path.isfile(fname)):
        loadNeeded = True
        break
    fname = "medianpressure_%s.txt" % months[m]
    if (not os.path.isfile(fname)):
        loadNeeded = True
        break
    fname = "medianhumidity_%s.txt" % months[m]
    if (not os.path.isfile(fname)):
        loadNeeded = True
        break

calculationNeeded = False
for m in xrange(0, len(months)):
    for f in xrange(0, len(frequencies)):
        fname = "opacity_%s_%.0f.txt" % (months[m], frequencies[f])
        if (not os.path.isfile(fname)):
            calculationNeeded = True
            break
        fname = "atmostemp_%s_%.0f.txt" % (months[m], frequencies[f])
        if (not os.path.isfile(fname)):
            calculationNeeded = True
            break
        fname = "atmostrans_%s_%.0f.txt" % (months[m], frequencies[f])
        if (not os.path.isfile(fname)):
            calculationNeeded = True
            break

# Read in the data.
if (loadNeeded):
    print "reading temperature"
    temperature = {}
    if (os.path.isfile(".site.environment.weather.Temperature.txt")):
        d = np.loadtxt(".site.environment.weather.Temperature.txt")
        temperature['time'] = []
        temperature['value'] = []
        for r in d:
            if (r[1] < 50):
                temperature['time'].append(datetime.datetime.fromtimestamp(int(r[0])))
                temperature['value'].append(r[1])

    print "reading pressure"
    pressure = {}
    if (os.path.isfile(".site.environment.weather.Pressure.txt")):
        d = np.loadtxt(".site.environment.weather.Pressure.txt")
        pressure['time'] = []
        pressure['value'] = []
        for r in d:
            pressure['time'].append(datetime.datetime.fromtimestamp(int(r[0])))
            pressure['value'].append(r[1] * 100.0)

    print "reading humidity"
    humidity = {}
    if (os.path.isfile(".site.environment.weather.RelHumidity.txt")):
        d = np.loadtxt(".site.environment.weather.RelHumidity.txt")
        humidity['time'] = []
        humidity['value'] = []
        for r in d:
            humidity['time'].append(datetime.datetime.fromtimestamp(int(r[0])))
            humidity['value'].append(r[1])

data = {
    'temperature': {},
    'pressure': {},
    'humidity': {}
    }
medianData = {
    'temperature': {},
    'pressure': {},
    'humidity': {},
    'opacity': {},
    'atmosphericTemperature': {},
    'atmosphericTransmissivity': {}
    }
medianInterval = 30.0 / 1440.0 # 10 minutes

if (loadNeeded):
    oldDate = datetime.datetime(startYear, 1, 1, 0, 0, 0)
    oldi = 0
    oldj = 0
    oldk = 0
    for y in xrange(startYear, endYear + 2):
        data['temperature'][str(y)] = {}
        data['pressure'][str(y)] = {}
        data['humidity'][str(y)] = {}
        print "= %d =" % y
        for m in xrange(0, len(months)):
            if (y == 2010 and m == 0):
                continue
            checkDate = datetime.datetime(y, m + 1, 1, 0, 0, 0)
            monthName = months[(m - 1)]
            ty = str(oldDate.year)
            data['temperature'][ty][monthName] = {
                'time': [], 'value': [], 'offset': [] }
            for i in xrange(oldi, len(temperature['time'])):
                if (temperature['time'][i] < checkDate):
                    data['temperature'][ty][monthName]['time'].append(temperature['time'][i])
                    data['temperature'][ty][monthName]['offset'].append(((temperature['time'][i] - oldDate).total_seconds()) / 86400.0)
                    data['temperature'][ty][monthName]['value'].append(temperature['value'][i])
                else:
                    oldi = i
                    break
            
            data['pressure'][ty][monthName] = {
                'time': [], 'value': [], 'offset': [] }
            for j in xrange(oldj, len(pressure['time'])):
                if (pressure['time'][j] < checkDate):
                    data['pressure'][ty][monthName]['time'].append(pressure['time'][j])
                    data['pressure'][ty][monthName]['offset'].append(((pressure['time'][j] - oldDate).total_seconds()) / 86400.0)
                    data['pressure'][ty][monthName]['value'].append(pressure['value'][j])
                else:
                    oldj = j
                    break
                
            data['humidity'][ty][monthName] = {
                'time': [], 'value': [], 'offset': [] }
            for k in xrange(oldk, len(humidity['time'])):
                if (humidity['time'][k] < checkDate):
                    data['humidity'][ty][monthName]['time'].append(humidity['time'][k])
                    data['humidity'][ty][monthName]['offset'].append(((humidity['time'][k] - oldDate).total_seconds()) / 86400.0)
                    data['humidity'][ty][monthName]['value'].append(humidity['value'][k])
                else:
                    oldk = k
                    break
            oldDate = checkDate

for m in xrange(0, len(months)):
    monthName = months[m]
    print "+ %s +" % monthName
    medianData['temperature'][monthName] = {
        'value': [], 'offset': [] }
    medianData['pressure'][monthName] = {
        'value': [], 'offset': [] }
    medianData['humidity'][monthName] = {
        'value': [], 'offset': [] }
    if (loadNeeded == False):
        fname = "mediantemperature_%s.txt" % monthName
        print "Loading %s" % fname
        mt = np.loadtxt(fname)
        medianData['temperature'][monthName]['offset'] = [ r[0] for r in mt ]
        medianData['temperature'][monthName]['value'] = [ r[1] for r in mt ]
        fname = "medianpressure_%s.txt" % monthName
        print "Loading %s" % fname
        mt = np.loadtxt(fname)
        medianData['pressure'][monthName]['offset'] = [ r[0] for r in mt ]
        medianData['pressure'][monthName]['value'] = [ r[1] for r in mt ]
        fname = "medianhumidity_%s.txt" % monthName
        print "Loading %s" % fname
        mt = np.loadtxt(fname)
        medianData['humidity'][monthName]['offset'] = [ r[0] for r in mt ]
        medianData['humidity'][monthName]['value'] = [ r[1] for r in mt ]
    
    medianData['opacity'][monthName] = {}
    medianData['atmosphericTemperature'][monthName] = {}
    medianData['atmosphericTransmissivity'][monthName] = {}
    for f in xrange(0, len(frequencies)):
        fstr = "f" + str(int(frequencies[f]))
        medianData['opacity'][monthName][fstr] = {
            'value': [], 'offset': [] }
        medianData['atmosphericTemperature'][monthName][fstr] = {
            'value': [], 'offset': [] }
        medianData['atmosphericTransmissivity'][monthName][fstr] = {
            'value': [], 'offset': [] }
    if (loadNeeded):
        print "Calculating medians"
        oldis = []
        oldjs = []
        oldks = []
        for y in xrange(startYear, endYear + 1):
            oldis.append(0)
            oldjs.append(0)
            oldks.append(0)
        for o in np.arange(0, ndays[m], medianInterval):
            nxo = o + medianInterval
            ilist = []
            jlist = []
            klist = []
            for y in xrange(startYear, endYear + 1):
                #print "searching temp from %d to %d" % (oldis[(y - startYear)], len(data['temperature'][str(y)][monthName]['offset']))
                for i in xrange(oldis[(y - startYear)], len(data['temperature'][str(y)][monthName]['offset'])):
                    if (data['temperature'][str(y)][monthName]['offset'][i] < nxo):
                        ilist.append(data['temperature'][str(y)][monthName]['value'][i])
                        oldis[(y - startYear)] = i + 1
                    else:
                        break
                for j in xrange(oldjs[(y - startYear)], len(data['pressure'][str(y)][monthName]['offset'])):
                    if (data['pressure'][str(y)][monthName]['offset'][j] < nxo):
                        jlist.append(data['pressure'][str(y)][monthName]['value'][j])
                        oldjs[(y - startYear)] = j + 1
                    else:
                        break
                for k in xrange(oldks[(y - startYear)], len(data['humidity'][str(y)][monthName]['offset'])):
                    if (data['humidity'][str(y)][monthName]['offset'][k] < nxo):
                        klist.append(data['humidity'][str(y)][monthName]['value'][k])
                        oldks[(y - startYear)] = k + 1
                    else:
                        break
            medianData['temperature'][monthName]['offset'].append(o + medianInterval / 2.0)
            medianData['temperature'][monthName]['value'].append(np.median(np.array(ilist)))
            medianData['pressure'][monthName]['offset'].append(o + medianInterval / 2.0)
            medianData['pressure'][monthName]['value'].append(np.median(np.array(jlist)))
            medianData['humidity'][monthName]['offset'].append(o + medianInterval / 2.0)
            medianData['humidity'][monthName]['value'].append(np.median(np.array(klist)))

    if (calculationNeeded):
        print "calculating atmospheric properties"
        for i in xrange(0, len(medianData['temperature'][monthName]['value'])):
            r = calcOpacity(hfrequencies, math.radians(90.0),
                            (medianData['temperature'][monthName]['value'][i] + 273.15),
                            medianData['pressure'][monthName]['value'][i],
                            (medianData['humidity'][monthName]['value'][i] / 100.0))
            for f in xrange(0, len(frequencies)):
                fstr = "f" + str(int(frequencies[f]))
                medianData['opacity'][monthName][fstr]['offset'].append(medianData['temperature'][monthName]['offset'][i])
                medianData['atmosphericTemperature'][monthName][fstr]['offset'].append(medianData['temperature'][monthName]['offset'][i])
                medianData['atmosphericTransmissivity'][monthName][fstr]['offset'].append(medianData['temperature'][monthName]['offset'][i])
                medianData['opacity'][monthName][fstr]['value'].append(r['tau'][f])
                medianData['atmosphericTemperature'][monthName][fstr]['value'].append(r['Tb'][f])
                medianData['atmosphericTransmissivity'][monthName][fstr]['value'].append(r['fac'][f])
        
if (loadNeeded):
    fig = plt.figure(figsize=(8.267, 11.692))
    for i in xrange(0, len(months)):
        plt.subplot(6, 2, (i + 1))
        plt.xlabel('Days since month start')
        plt.ylabel('Temperature (C)')
        plt.title(months[i])
        for j in xrange(startYear, endYear + 1):
            plt.plot(data['temperature'][str(j)][months[i]]['offset'],
                     data['temperature'][str(j)][months[i]]['value'], linecols[(j - startYear)], label=str(j))
        plt.axis([0, 31, -5, 50])
    plt.tight_layout()
    plt.savefig('temperatureplot.pdf')

fig = plt.figure(figsize=(8.267, 11.692))
for i in xrange(0, len(months)):
    plt.subplot(6, 2, (i + 1))
    plt.xlabel('Days since month start')
    plt.ylabel('Temperature (C)')
    plt.title(months[i])
    plt.plot(medianData['temperature'][months[i]]['offset'],
             medianData['temperature'][months[i]]['value'], linecols[0])
    fname = 'mediantemperature_%s.txt' % months[i]
    np.savetxt(fname, (medianData['temperature'][months[i]]['offset'],
                       medianData['temperature'][months[i]]['value']))
    plt.axis([0, 31, -5, 50])
plt.tight_layout()
plt.savefig('mediantemperatureplot.pdf')

if (loadNeeded):
    fig = plt.figure(figsize=(8.267, 11.692))
    for i in xrange(0, len(months)):
        plt.subplot(6, 2, (i + 1))
        plt.xlabel('Days since month start')
        plt.ylabel('Pressure (Pa)')
        plt.title(months[i])
        for j in xrange(startYear, endYear + 1):
            plt.plot(data['pressure'][str(j)][months[i]]['offset'],
                     data['pressure'][str(j)][months[i]]['value'], linecols[(j - startYear)], label=str(j))
        plt.axis([0, 31, 95000, 105000])
    plt.tight_layout()
    plt.savefig('pressureplot.pdf')

fig = plt.figure(figsize=(8.267, 11.692))
for i in xrange(0, len(months)):
    plt.subplot(6, 2, (i + 1))
    plt.xlabel('Days since month start')
    plt.ylabel('Pressure (Pa)')
    plt.title(months[i])
    plt.plot(medianData['pressure'][months[i]]['offset'],
             medianData['pressure'][months[i]]['value'], linecols[0])
    fname = 'medianpressure_%s.txt' % months[i]
    np.savetxt(fname, (medianData['pressure'][months[i]]['offset'],
                       medianData['pressure'][months[i]]['value']))
    plt.axis([0, 31, 95000, 105000])
plt.tight_layout()
plt.savefig('medianpressureplot.pdf')

if (loadNeeded):
    fig = plt.figure(figsize=(8.267, 11.692))
    for i in xrange(0, len(months)):
        plt.subplot(6, 2, (i + 1))
        plt.xlabel('Days since month start')
        plt.ylabel('Humidity (%)')
        plt.title(months[i])
        for j in xrange(startYear, endYear + 1):
            plt.plot(data['humidity'][str(j)][months[i]]['offset'],
                     data['humidity'][str(j)][months[i]]['value'], linecols[(j - startYear)], label=str(j))
        plt.axis([0, 31, 0, 100])
    plt.tight_layout()
    plt.savefig('humidityplot.pdf')

fig = plt.figure(figsize=(8.267, 11.692))
for i in xrange(0, len(months)):
    plt.subplot(6, 2, (i + 1))
    plt.xlabel('Days since month start')
    plt.ylabel('Humidity (%)')
    plt.title(months[i])
    plt.plot(medianData['humidity'][months[i]]['offset'],
             medianData['humidity'][months[i]]['value'], linecols[0])
    fname = 'medianhumidity_%s.txt' % months[i]
    np.savetxt(fname, (medianData['humidity'][months[i]]['offset'],
                       medianData['humidity'][months[i]]['value']))
    plt.axis([0, 31, 0, 100])
plt.tight_layout()
plt.savefig('medianhumidityplot.pdf')

fig = plt.figure(figsize=(8.267, 11.692))
for i in xrange(0, len(months)):
    plt.subplot(6, 2, (i + 1))
    plt.xlabel('Days since month start')
    plt.ylabel('Opacity')
    plt.title(months[i])
    for j in xrange(0, len(frequencies)):
        fstr = "f" + str(int(frequencies[j]))
        plt.plot(medianData['opacity'][months[i]][fstr]['offset'],
                 medianData['opacity'][months[i]][fstr]['value'], linecols[j])
        fname = 'opacity_%s_%.0f.txt' % (months[i], frequencies[j])
        np.savetxt(fname, (medianData['opacity'][months[i]][fstr]['offset'],
                           medianData['opacity'][months[i]][fstr]['value']))
    #plt.axis([0, 31, 0, 100])
plt.tight_layout()
plt.savefig('opacityplot.pdf')

fig = plt.figure(figsize=(8.267, 11.692))
for i in xrange(0, len(months)):
    plt.subplot(6, 2, (i + 1))
    plt.xlabel('Days since month start')
    plt.ylabel('Atmospheric Temperature')
    plt.title(months[i])
    for j in xrange(0, len(frequencies)):
        fstr = "f" + str(int(frequencies[j]))
        plt.plot(medianData['atmosphericTemperature'][months[i]][fstr]['offset'],
                 medianData['atmosphericTemperature'][months[i]][fstr]['value'], linecols[j])
        fname = 'atmostemp_%s_%.0f.txt' % (months[i], frequencies[j])
        np.savetxt(fname, (medianData['atmosphericTemperature'][months[i]][fstr]['offset'],
                           medianData['atmosphericTemperature'][months[i]][fstr]['value']))
    #plt.axis([0, 31, 0, 100])
plt.tight_layout()
plt.savefig('atmospherictemperatureplot.pdf')

fig = plt.figure(figsize=(8.267, 11.692))
for i in xrange(0, len(months)):
    plt.subplot(6, 2, (i + 1))
    plt.xlabel('Days since month start')
    plt.ylabel('Atmospheric Temperature')
    plt.title(months[i])
    for j in xrange(0, len(frequencies)):
        fstr = "f" + str(int(frequencies[j]))
        plt.plot(medianData['atmosphericTransmissivity'][months[i]][fstr]['offset'],
                 medianData['atmosphericTransmissivity'][months[i]][fstr]['value'], linecols[j])
        fname = 'atmostrans_%s_%.0f.txt' % (months[i], frequencies[j])
        np.savetxt(fname, (medianData['atmosphericTransmissivity'][months[i]][fstr]['offset'],
                           medianData['atmosphericTransmissivity'][months[i]][fstr]['value']))
    #plt.axis([0, 31, 0, 100])
plt.tight_layout()
plt.savefig('atmospherictransmissivityplot.pdf')
